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Abstract 



The material in this article was presented in five hours of lectures to the 2011 
Tenerife Winter School. The School's theme was "Secular Evolution of Galaxies" 
and my task was to present the underlying stellar-dynamical theory. Other 
lecturers were speaking on the role of bars and chemical evolution, so these 
topics are avoided here. The material starts with an account of the connections 
between isolating integrals, quasiperiodicity and angle-action variables - these 
variables played a prominent and unifying role throughout the lectures. This 
leads on to the phenomenon of resonant trapping and how this can lead to 
chaos in cuspy potentials and phase-space mixing in slowly evolving potentials. 
Surfaces of section and frequency analysis are introduced as diagnostics of phase- 
space structure. Real galactic potentials include a fluctuating part that drives 
the system towards unattainable thermal equilibrium. Two-body encounters are 
only one source of fluctuations, and all fluctuations will drive similar evolution. 
The orbit-averaged Fokker-Planck equation is derived, as are relations that hold 
between the second-order diffusion coefficients and both the power spectrum of 
the fluctuations and the flrst-order diffusion coefficients. From the observed 
heating of the solar neighbourhood we show that the second-order diffusion 
coefficients must scale as ~ J^/^. We show that periodic spiral structure shifts 
angular momentum outwards, heating at the Lindblad resonances and mixing 
at corotation. The equation that would yield the normal modes of a stellar disc 
is first derived and then used to discuss the propagation of tightly-wound spiral 
waves. The winding up of such waves is described and explains why cool stellar 
discs are responsive systems that amplify ambient noise. An explanation is 
offered of why the Lin-Shu-Kalnajs dispersion relation and even global normal- 
mode calculations provide a very incomplete understanding of the dynamics of 
stellar discs. 
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1.1 Orbits 

Even a gas-free galaxy is a formidably complex dynamical system: ~ 10^^ 
stars and inconceivably more WIMPS are strongly coupled by their mutual 
gravitational field. As in any branch of theoretical physics we make progress 
by simplifying. The key simplification is that the motion of all these particles 
would be very similar if calculated in the gravitational field of an imaginary 
smoothed-out version of the galaxy. That is, we smear the mass of each particle 
over something like an interstellar distance and compute the gravitational field 
from this smooth mass distribution. Given the usefulness of this approximate 
gravitational field, our first task becomes to understand how particles move in 
smooth gravitational fields of this kind. 
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Figure 1.1: The motion in the {R,z) plane of a particle that orbits in an ax- 
isymmetric gravitational potential. The orbits have the same values of E and 
Lz but are nevertheless very different because they have different values of the 
mysterious "third integral" . 



1.1.1 Quasiperiodicity and isolating integrals 

The equations of motion of a particle in a typical gravitational potential $ are 
readily integrated numerically. If the potential is axisymmetric, the component 
of angular momentum about the galaxy's symmetry axis, Lz, is conserved, and 
we can eliminate the azimuthal angle (j) and its time derivative from the equa- 
tions of motion, leaving us with two coupled ordinary differential equations to 
integrate: 

a<I>cff ., d^eS ^^-^^ 



R 



dR 



dz 



where 



<i>eff(i?,z) = <i>(i?,z) + ^. 



(1.2) 



Fig. 11.11 shows a couple of typical orbits obtained in this way. They have a 
nicely regular pattern of paths running to and fro. Given sufficient time such 
an orbit will carry the particle to any given point inside the bounding envelope, 
and when it reaches that point the particle will be moving with one of four 
velocities: R = and z = ±Vz, where vr and Vz are numbers that depend 
on the orbit and the location. 

If the time series R{t) or z{t) obtained in this way is Fourier transformed, one 
finds that the spectrum consists of discrete lines, and moreover, any frequency 
appearing in the spectrum can be expressed as an integer linear combination of 
two fundamental frequencies fi^ and ilz- That is. 



R{t) = ^ i?„ cos(n • nt + ipRn) 



z(t) 



Z„ cos(n • Jit 



-V'.n),(1.3) 



where 



Q. ^ {VLii,VLz) and n^{nii,nz) 



(1.4) 



with rifl and Uz being (possibly negative) integers. Thus the orbit is charac- 
terised by the two fundamental frequencies and countably many amplitudes 
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Rn,Zn and phases ipB.n,ipzn- An orbit with these properties is said to be 
quasiperiodic. 

As the particle moves in the (R, z) plane, it rotates about the symmetry 
axis: its azimuthal coordinate (p can be found by simple quadrature: 

^{t)^m+ J^dt^. (1.5) 

The time series of (j) will consist of sinusoidal contributions at discrete frequencies 
superimposed on a secularly increasing value, and the rate of the secular increase 
defines a third fundamental frequency $7^. Thus three fundamental frequencies 
will be required to express the frequency of any line in the spectrum of say x(t) 
as an integer linear combination of fundamental frequencies. 

The result of integrating the equations of motion in potentials that are not 
axisymmetric cannot be displayed as nicely as in Fig. 11.11 But when you Fourier 
transform the time series x{t) of the coordinates of such orbits, you usually 
find that the spectra consist of lines that can be expressed as integer linear 
combinations of three fundamental frequencies. So a large fraction of orbits in 
galaxy-like potentials are quasiperiodic. 

An integral of motion is a function /(x, v) on phase space that returns the 
same number no matter at what point along an orbit you evaluate it: 

/(x(t),v(t)) = constant. (1.6) 

We say that / is an isolating integral if the set of phase-space points that 
satisfy equation l|1.6p for some given value on the right side defines a smooth 
five-dimensional subset of phase space. For example, the energy i?(x, v) = 
iz)^ -|- <i>(x) is an isolating integral. 

If an orbit is quasiperiodic, one may show that it admits at least three func- 
tionally independent isolating integrals. We can take E to be one of these, and 
in the case of an axisymmetric potential, — xvy — yvx can be taken to be an- 
other of them. Only in exceptional potentials do we have an analytic expression 
for a third isolating integral, bu t its exi s tence is assured by the quasiperiodic 



nature of orbits - for a proof see lArnoldl (|1978[ ) 



1.1.2 Angles and actions 

It often happens that the key to solving a physics problem is to identify the 
coordinate system that's best suited to that problem. Hamiltonian mechanics 
allows us to use an extremely wide range of coordinates for phase space with 
ease - all "canonical coordinates" are intrinsically equal. In the following we 
shall denote by (x, p) a canonical system made up of coordinate Xi that gives 
the star's position and its canonically conjugate momentum p^. In the simplest 
case Xi IS a Cartesian coordinate and pi = Xi is its rate of change, but in 
some instances pi ^ Xi. Given that three isolating integrals li exist, a shrewd 
question to ask is "is there a canonical coordinate system in which the 7^, or 
some functions of them, act as momenta?" Since any function J{Ii, I2, 13) of 
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Figure 1.2: Trajectories generated by an action (left) and by some other integral 
(right). In the left panel each horizontal curve joins points on the left and right 
boundaries that have been identified, whereas in the right panel the leftmost and 
rightmost points of these curves are at logically distinct points. Consequently, 
in the right panel the periodic continuation of each horizontal curve would not 
overlie the portion plotted. The vertical curves exhibit the same phenomenon. 
Consequently these curves define a global coordinate system only in the left 
panel. 



the li will itself be an integral of motion, we usefully increase our chances of 
our enquiry having the answer "yes" by opening the enquiry up to functions of 
our original integrals. 

It turns out that only a highly restricted group of integrals are capable of 
playing the role of momenta. We call such integrals actions and reserve for 
them the letter J, so Jj(x, p) is the jth action of an orbit. 

To understand why only special integrals can act as momenta, we should 
consider the subset M of points in phase space that satisfy 

/i(x, p) = constant for i — 1,2,3. (1-7) 

These equations impose three restrictions on the six phase-space coordinates 
(x, p). So the set M is a three-dimensional subset of phase space to which 
the particle is confined for all time. M is a subset of the energy hypersurface 
i?(x, p) — E, where H is the Hamiltonian function and E is the particle's 
energy. If the orbit is bound, as we shall assume, the energy hypersurface is 
compact, so M is compact also. In practice M will also be a connected set, and 
it can then be shown that it is diffeomorphic to a 3-torus (see §§49,50 of lArnoldl . 



19781 for a proof). 

What is a 3-torus? Think of a room with each point on the fioor identified 
with the point on the ceiling that is vertically above it, each point on the left 
wall identified with the corresponding point on the right wall, and with points 
on the front and back walls similarly identified. Position within this room (M) 
is specified by the values of the coordinates 6i, 62 and 6*3 that are canonically 
conjugate to Ji, J2 and J3, whose values specify which room (torus) we are in. 
They do so in a way that can be understood geometrically. 
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In Hamiltonian mechanics a fundamental role is played by the Poincare 
invariant P. This is a number that we assign to any two-surface S in phase 
space through 

P{S) = J2[ dxidpi, (1.8) 
■ Js 

so P{S) is the sum of the areas of the projections of S onto all the planes formed 
by each coordinate and its conjugate momentum pj. It's called an invariant 

because you get the same mimbcr no matter what canonical coordinates you 
use. Consequently, given any two-surface S & M we can replace {xi,Pi) by 
{9i, Ji) and write 

P{S) = Y. I AOidJi (1.9) 



s 

and this vanishes because every point of S has the same value of Jj. In view of 
this fact we say M is null. It follows from the nullness of M that the value of 
any line integral J dx • p through M between two given endpoints is the same 
for any path between those points that can be continuously distorted into some 
standard path without leaving M. To see why, take the difference between the 
integrals along two different paths Fi and r2 between a given pair of points. 
This difference of integrals is identical with the line integral along the closed 
path Fi — F2 in which we go out on Fi and back on F2. By Green's theorem 
this closed line integral is equal to the Poincare invariant of the 2-surface that 
r2 sweeps out as it is deformed into Fi. But this Poincare invariant vanishes, 
so the original line integrals were equal. Consider now the line integral / dx • p 
along the path from the front wall to the corresponding point on the back wall. 
We choose to evaluate this integral using the {6, J) coordinates, and to take the 
path on which 02 = 63 = constant. Then we have 



j dx-p = j dej = Ji j 



dOi (1.10) 



so Ji is equal to / dx • p divided by the amount by which 9i increments as 
we cross the room. We choose to scale the actions so 9i increments by 27r on 
crossing the room, so 

= i /r/^ • P = i ^ '^''^P^' ^^-^^^ 

where F^ is any path that takes us once across the room. 

If we want to find the values of the ordinary phase-space coordinates (x, p) 
at the point 9 we need to do the following: (i) Choose a point (xo,po) in the 
room and declare it to be = 0. Now integrate from the initial conditions 
x = Xo, p = po when 9i = the coupled o.d.e.s 

^ = [^i,^i], = (1-12) 
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where [..] denotes a Poisson bracket, to discover the (x, p) coordinates of the 
points 9 = (01, 0, 0). (ii) Starting from any of these points we can integrate the 
o.d.e.s 

^ = [^.,^2], ^ = (1.13) 

to discover the (x, p) coordinates of the points 6 = {61,62,0). (iii) Starting 
from any of the last-mentioned points we can integrate the third set of o.d.e.s 
to find the (x, p) coordinates of a general point 6. The key property of actions 
is that if we integrate say the first set of o.d.e.s from a point 9a on one wall 
to the point 9b at which the trajectory hits the opposite wall, we find that 
9b = 9a + (27r, 0, 0) (Fig. 11.21 left). That is the trajectories generated by actions 
carry you from a point on one wall to the point on the opposite wall that we 
have identified with our starting point. If we form some other integrals /i(J) by 
adopting three non-trivial functions of the actions, the trajectories generated 
by the Jj will not carry us from a point on a wall to the point on the opposite 
wall that has been identified with it (Fig. ll.2l right). Consequently, the variables 
that are canonically conjugate to the Ii cannot form a global coordinate system 
(although they can be used as coordinates locally). 

In remembrance of the fact that 9a and 9b = 9a + (27r, 0,0) correspond to 
the same point in phase space, the canonically conjugate coordinates of actions 
are called angle variables. To see what flexibility we have in the choice of 
actions and angles, we observe that if we have new variables 9' that are related 
to our original angle variables by linear equations 

6> = N-6>', (1.14) 

where N is a matrix with integer entries, then if any of the 6[ is incremented by 
27r, all the 6j will change by 2m7r, so we will return to the same place in phase 
space. Hence the 9' provide a global coordinate system for a torus as effectively 
as the 9. To discover what actions correspond to the 9' , we write down the 
generating function of the canonical transformation {9, J) -r-> {9', J') 

5(0',J) J,7V„-0;-. (1.15) 
It is straightforward to check that S generates equation (11.141) and we have also 

'^i = a^ = E^^^^- (1-16) 

Thus the new actions are integer linear combinations of the old actions. 

Any function on phase space can be expressed as a function of the angles and 
actions {9, J). In particular the Hamiltonian can be expressed in this form. But 
while the angle variables vary as a particle orbits, its energy does not. So the 
Hamiltonian cannot depend on 9. Hence the Hamiltonian is a function H(3) of 
the actions only. This fact makes the equations of motion of the angle variables 
trivial: 

dH ^ . . 
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Figure 1.3: A particle moves tlirough its room (torus) on a straight line that 
repeatedly hits a wall and reappears at the corresponding point on the opposite 
wall. If the components of the slope vector are rationally related, the line even- 
tually closes on itself (right). In general it does not close and comes arbitrarily 
close to every point in the room (left). Black dots mark the start and end of 
the plotted portion of the trajectory. 

Since H depends on J only, so must its partial derivatives fli. Since J is constant 
along the orbit, il; is too, so we can immediately integrate equation (jl.l7l) : 

9,{t) ^ e,{o) + n,t. (1.18) 

Thus the particle moves through its room at constant speed along straight lines 
whose slope is given by S7 = dH/dJ. 

Usually the frequencies are incommensurable and then the line will even- 
tually come arbitrarily close to every point in the room (Fis. 11.31 left). 

Since the ordinary phase-space coordinates are periodic functions of the angle 
variables with period 2tt they can be expanded in Fourier series 

a:;,(6/) = ^X,„cos(n-6/ + ?A,n), (1.19) 

n 

where the sum is over vectors with integer components and the Xin and ipia are 
constant amplitudes and phases. When we substitute into this expression our 
solution (|1.18p of the equations of motion, we obtain 

a:^(t) =^X,nCOs(n-rJi + V'-n), where = V'in + n • 6/,(0). (1.20) 

n 

Equation (|1.3p from which we started is an instance of this equation. Thus we 
have come full circle from the empirical fact that numerically integrated orbits 
are quasiperiodic to the fact that such spectra are a necessary consequence of 
these orbits having three isolating integrals. 

In the generic case of incommensurable frequencies, we have a useful result, 
the time-averages theorem: when the frequencies are incommensurable, the 
fraction of its time a particle spends in a subset V of the torus is Jy d'^0/(27r)'^. 
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From this the strong Jeans theorem follows: in a steady- state galaxy, the 
density of stars is uniform within incommensurable tori, so the density of stars 
in phase space is a function /(J) of the actions only. We call / the distribution 
function, often abbreviated DF. 

It is useful to think of orbits as points in the three-dimensional space, action 
space, that has the three actions as is Cartesian coordinates. The strong Jeans 
theorem tells us that galactic equilibria are simply distributions of stars in this 
easily imagined space. 



1.1.3 Epicycle approximation 

It's instructive to examine a very useful model at this point. A star in an 
axisymmetric potential moves in the effective potential (jl.2p . It is physically 
obvious that the minimum of ^cS occurs at (i?g, 0), where the guiding-centre 
radius i?g is the radius of the circular orbit with the given angular momentum 
Lz- Since this is the potential's minimum, a Taylor expansion of $off around it 
will contain no linear terms and be of the form 

*eff(i?g + x, z) = $cff(i?g, 0) + ^K^x' + Iv'z' + • • • (1.21) 

Stars on sufficiently circular orbits will be confined to the region in which we 
need retain only the first three terms in this series. Consequently, their radial 
and vertical motions will be harmonic. The frequencies of these oscillations 
are il^ = k, the epicycle frequency and = the vertical epicycle 
frequency. The solution to the equations of motion is 

i?(i) = i?g + Xcos(Ki + ?/'r) ; z{t) = Zco's.{vt + i}z), (1.22) 

where X, V'r, ■Z^ and -02 are all constants. Clearly we can set 

6lr = Kt + Vr ; 9z = vt^'\l}z. (1.23) 

Then calculating p^j — R and evaluating ^ dRpn we can show that X = 
y^2Jr/ K, and similarly that Z =^ ■y/2 J^/z^, so in the epicycle approximation 
the connection between ordinary coordinates and angle-action variables is 



2J 2J 
+ J -JLcosOr ; z^d—j-cosOz. (1.24) 



1.2 Resonances 

The fundamental frequencies Jl are generally non-trivial functions of the actions, 
so for certain tori a resonance condition N • 51 = will apply. When we 
multiply the resonance condition by t and use equation (ll.lSp . we obtain 

N • e{t) = constant. (1.25) 

This equation inspires us to make a canonical transformation to new angles and 
actions (0', J') such that N • 6 becomes one of the new angles, say O'l- It does 
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not much matter what we adopt for 62 and d'^\ 0'2 = O2 and ^3 = ^3 works fine. 
Then our generating function becomes 

S'(0,J') = J(N-6/+ J^6i2 + J36I3, (1.26) 

so the new angles are 

9[^gjr^N-e, 9'2^—^^e,, 9',^e,. (1.27) 

Since 9'i does not evolve in time, the star explores only a two-dimensional set 
of its three-dimensional torus. While a star on a resonant torus does not come 
arbitrarily close to every point on its torus, the bigger the numbers Nj are, the 
more effectively it samples its torus, and the less likely it is that the resonance 
condition will be dynamically important. 

1.2.1 Perturbation theory 

The importance of resonances only emerges when we consider the effects of 
adding a small perturbing Hamiltonian h{6,J) to the Hamiltonian Ho{3) for 
which we have found angle-action variables (0, J). We use these coordinates to 
study the motion of a particle in the full Hamiltonian H = Hq + h. Hamilton's 
equations read 

where JIq — dHo/dJ. The perturbing Hamiltonian, like any function on phase 
space, is a periodic function of the angles, so we can Fourier expand it: 

/i(6>,J) =^/i„(J)cos(n-6>-|-V'„). (1.29) 
11 

With this expansion, the equation of motion of J is 

j = ^ n/i„(J) sin(n • + V^„) = ^ n/i„(J) sin(n -nt + iP'J. (1.30) 



So long as n • n 7^ 0, the time-averaged value of J vanishes and we expect J 
simply to oscillate slightly around its unperturbed value. But if a resonance 
condition is nearly satisfied, N • $7 ~ 0, the argument of one or more of the sines 
will change very slowly, and the cumulative change in J can be significant even 
for very small /in- Thus resonances permit small forces to act in the same sense 
for long times, and hence cause qualitative changes in behaviour. This is one of 
the fundamental principles of physics. 

Mathematically, we use the new angle variables 6' defined by equation (|1.27p 
and their conjugate actions 

Jl^^^J'l^l, J2 = — - J(7V2 + J2, J3=J[N3 + 4. (1.31) 
uui OU2 
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Next we express /i as a Fourier series in the new angle variables and discard 
all terms that involve 62 or 9'^ on the grounds that they oscillate too rapidly to 
have a significant impact on the dynamics. Since our approximated Hamiltonian 
depends on neither 62 nor ^3, the conjugate actions J2 and J3 will be constants 
of motion. The only non-trivial equations of motion are now 

j'l^T. -^^'n^J'^) + ^(«,o,o)), (1-32) 

n 

where 

BH ( T'l B J 

" BJ' " ^ ^"'W " = N • f^o (1.33) 

and we have suppressed the dependence of Qqi and h'^^ on J2 and J3 because 
the latter are mere constants. We have reduced the particle's motion in six- 
dimensional phase space to motion in the (O'l, J() plane. 

We differentiate with respect to time the first of equations (|1.32p 

BO' /B'^h' Bh' \ 

(1-34) 

We can neglect the sum in this equation because each of its terms is the product 
of a derivative of /i^, which is small, and either J{, which is of the same order, 
or ^5.1 which is also small because il'^i vanishes at the resonance. Therefore we 
can dramatically simplify the 9'i equation of motion to 

BQ' BO' 

If we approximate dO'Q^/BJ[ and h'^ by their values on resonance, and retain 
only the term for n = 1 in the sum, we are left with the equation of motion of 
a pendulum 

= ~'Uj^sm0, (1.36) 

where 

BO' 

e = e[+i;[ and cj^ = -—f-h[. (1.37) 

(yu ^ 

A pendulum can move in two ways: at low energy its motion is oscillatory 
at an angular frequency that falls from uj at the lowest energies to zero at the 
critical energy above which the pendulum circulates. Consequently, equation 
p. 351) predicts that close to the resonance ("low energy") 9'i will oscillate. This 
is the regime of resonant trapping in which the particle librates around the 
underlying resonant orbit. At some critical distance from the resonance ("high 
energy" ) 6'i will start to circulate. Near the threshold energy, the rate at which 
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Figure 1.4: The phase plane of a pendulum. Curves of constant energy E 
(eq. 11.381) are plotted. The particle moves on these, from left to right in the 
upper half of the figure and from right to left in the lower half. 

6'i advances in time will be highly non-uniform, just as a pendulum that has 
only just enough energy to get over top dead centre slows markedly as it does 
so. As we move further and further from the resonance, the rate of advance of 
0'i becomes more and more uniform, and we gradually recover the unperturbed 
motion, in which the rate of advance of 6'i is strictly uniform. 

We can obtain a useful pictorial representation of this behaviour by deriving 
the energy invariant of equation ()1.36p . We multiply both sides by 6 to obtain 
an equation which states that 

E=\d'^-iJcose (1.38) 

is constant. Consequently, the particle moves in the {9, 0) plane along curves of 
constant E. These curves are shown in Fig. 11.41 The round contours near the 
centre of the figure show the motion of a particle that has been trapped by the 
resonance, while the contours at the top and bottom of the figure that run from 
9 = —TT to 6 = TT show the motion of a particle that continues to circulate. 

To appreciate the significance of resonances for galactic dynamics we must 
introduce the concept of a surface of section. This is a phase plane such as the 
(R,piij plane for motion in an axisymmetric potential. Each time the particle 
passes through the plane z = moving upward {pz > 0) we mark the surface of 
section with a dot at the current values of R and pr. If the orbit is quasiperiodic, 
the dots eventually trace out a curve. This curve is the intersection of the orbit's 
torus with the {R,pji) plane, and in fact the area inside it, J dRdpu is 27r times 
the orbit's radial action Jr. Fig. 11.51 shows surfaces of section for motion in 
two flattened axisymmetric potentials. In a given panel each curve is for a 
different orbit with the same energy and value of L^- Most of the curves move 
around a central point. The point itself is made by the shell orbit, which is 
two-dimensional and has = 0; in real space this orbit is a thin cylindrical 
shell that has a larger diameter at z = than at its top or bottom edges. Each 
curve around this point in Fig. ll.5l is generated by a three-dimensional orbit that 
forms an annulus of finite thickness. Fig. 11.11 shows cross sections through two 
such annuli. In Fig. 11.51 the longer an orbit's curve is, the thicker the annulus 
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Figure 1.5: Surfaces of section for motion in flattened isochrone potentials: the 
left panel is for the case of a mass distribution that has axis ratio q — 0.7, 
while the right panel is for q — 0.4. In the right panel we see resonant islands 
generated by the 1 : 1 resonance between the radial and vertical oscillations. 
No such island is evident in the left panel. 



and the smaller its vertical extent. The outermost curve in Fig. ll.Sl is generated 
by the orbit Jz = 0, which is confined to the plane z = 0. So in Fig. 11.11 the 
orbit in the left panel generates a larger curve in the left panel of Fig. 11.51 than 
does the orbit of the right panel of Fig. 11.11 

The right panel in Fig. 11.51 is for motion in the potential of a flatter galaxy 
than the left panel, and in this panel not all curves loop around the central 
point. Two resonant islands have appeared, formed by orbits that have been 
trapped by the resonance 

At this point one should be asking "so what's the perturbation in this case?" 
It isn't evident that there is one; all we did was integrate orbits in a perfectly 
well defined potential. However, we can imagine that our Hamiltonian is the 
sum of a Hamiltonian Hq that would generate a surface of section in which 
all curves simply looped around the central point, and a smaller Hamiltonian 
h = H — Ho, and to ascribe to h the trapping of orbits by the 1 : 1 resonance. 



Kaasalainen &: Binnev ( 1994 ) describe an algorithm that can be used to generate 



Hq and h. 

The resonance just discussed is rather a tame one that probably does not 
play a big role in galactic dynamics (but see !jl.2.3p . Resonances certainly play 
a big role in the dynamics of bars. The dynamics of a bar is strongly affected 
by the bar's pattern speed, ftp, the angular velocity at which the figure of the 
bar rotates. The bars we see at the centres of more than half spiral galaxies, 
including our own, are rapidly rotating in the sense that in the rotating frame 
of the bar the dynamics of most stars depends strongly on the Coriolis force 
2J7p X p. In cosmological simulations of the clustering of dark matter, most 
dark halos are barred and their pattern speeds are so small that the Corio- 
lis force is unimportant for a significant fraction of their constituent particles. 
The gravitational potentials of galaxies are approximately isothermal, so let's 
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Figure 1.6: Surfaces of section for the potential (|1.39p with axis ratio q ~ 0.8. 
Left: for ro = 0.1; right: for ro = 0.01. 



investigate motion in the potential 

$ = ii;2 In j^r^ + + ^ . (1.39) 

Here rg is a core radius, within which the potential tends towards that of a 
harmonic oscillator, and vq would be the circular speed at x ^ ro in the (x, y) 
plane if the axis ratio q were unity. We shall assume that q = 0.8, so the 
potential is slightly elongated along the x axis and for now restrict ourselves to 
motion in the {x, y) plane so we can use the device of a surface of section. 

Fig. 11.61 shows {x,px) surfaces of section for two values of tq. In the case 
ro = 0.1 on the left we can identify the curves of two types of orbit: the two 
bulls-eyes comprise the curves of loop orbits such as that depicted in the top 
left panel of Fig. 11.71 These orbits have a well-defined sense of rotation about 
the z axis. The curves that surround both bulls-eyes are generated by box 
orbits such as that depicted in the bottom left panel of Fig. 11.71 Box orbits 
extend to the potential's centre and do not involve rotation. In the surface of 
section on the right of Fig. 11.61 we see several resonant islands in the domain of 
the box orbits. These islands are made up of orbits trapped by the resonances 
fly = 2nx: (outermost island) and 2il,y — 3flx (further in). The right panels of 
Fig. 11.71 show the appearance of a couple of these orbits in real space. As this 
experiment suggests, the more cuspy a triaxial mass distribution is, the larger is 
the fraction of pha se that is occupied by res onant box orbits, or boxlets as they 
have been called ( Merritt fc ValhiriL 1 19991) . Why resonant trapping becomes 
more important as the matter distribution becomes more cuspy is not clear. A 
fact that is probably relevant is that boxlets generally avoid the galactic centre, 
while boxes do not. Given that massive black holes reside in galaxy centres, 
the tendency of boxlets to avoid the centre clearly has important astrophysical 
consequences by keeping stars safe from being eaten by the monster there. 
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Figure 1.7: A loop and a box orbit (left); two boxlets (right). The cross in the 
lower right panel marks the potential's centre. 
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Figure 1.8: Frequency-ratio diagram for the potential ()1.39p with q = 0.9 and 
92 = 0.7. 
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Figure 1.9: A model Galaxy potential. The cuspiness of the equipotential con- 
tours clearly show the extent to which the potential has been flattened by the 
mass of the disc. 

1.2.2 Frequency analysis 

We have seen how surfaces of section give valuable insight into the structure 
of phase space. Unfortunately, it is impracticable to construct a surface of sec- 
tion unless the motion is in two dimensions - or effectively so in the case of an 
axisymmetric potential. Frequency analysis provides a useful way of determin- 
ing which resonances play an important role when the potential is triaxial and 
the motion of particles is inherently three-dimensional. The technique consists 
of numerically integrating orbits and then Fourier transforming the time de- 
pendence of suitable coordinates. From the resulting spectra one identifies the 
fundamental frequencies, and determines two frequency ratios, such a fi^/Oj, 
and riz/Qy. Each orbit then generates a point in a frequency-ratio diagram 
such as that of Fig. 11.81 In this diagram the points of non-resonant orbits fall 
along a curvilinear grid, which reflects the systematic way in which the initial 
conditions explored phase space. The points of resonant orbits lie along straight 
lines. In the vicinity of these lines there is a deficit of points, because orbits in 
the underlying integrable potential that would have the nearly commensurable 
frequencies that correspond to these locations have been resonantly trapped so 
their principal frequencies exactly satisfy a resonance condition. Actually these 
orbits still have three independent fundamental frequencies, but one of these 
frequencies is the frequency of libration around the trapping orbit, and this fre- 
quency has to be dug out of the coordinates' spectra and is missed by the simple 
algorithm from used to identify the frequencies from which the frequency ratios 
were calculated. 
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1.2.3 Adiabatic deformation 

The Milky Way's disc has accumulated rather gradually over most of the Hubble 
time. As it grew, the Galaxy's potential must have deformed from being the 
nearly spherical potential of the dark halo, to a potential that is significantly 
pinched towards the plane (Fig. 11.91) . Individual orbits will have distorted in 
response to the distortion of the potential, but so long as they remained non- 
resonant, their actions will have been invariant: actions are adiabtatic invari- 
ants. This fact greatly facilitates the process of determining the response of a 
stellar system to adiabatic distortion of its potential because all we have to do 
is to move each star from its orbit in the original potential to the orbit with 
the same actions in the distorted potential. In particular, the structure of the 
distorted system depends only on the initial and final configurations, and not 
on which configurations it passed through in between. 

The story is more complex and interesting if resonant trapping is possible. 
When the Galaxy's potential was nearly spherical, every orbit had a value of 
fir that was bigger than that of either of its other two frequencies. Stars whose 
orbits are now confined within a couple of kiloparsecs of the equatorial plane 
have fir < fiz, the frequency associated with motion perpendicular to the plane. 
So these stars have at some point satisfied the resonance condition fi^ = fiz . In 
a flatted potential fi^/fiz is smallest for orbits that are confined to the equatorial 
plane. Therefore the resonance condition was first satisfied by these orbits. 

Fig. 11.51 shows surfaces of section for motion in a potential before and after 
the resonance condition fi^ = fiz is first satisfied: the islands visible in the right 
panel are made up of orbits trapped by this resonance. Note that the areas of 
the curves in the left panel are 27rJr, so they do not change as the potential 
flattens. 

The resonance condition is first satisfied by the orbit that is confined to the 
equatorial plane; in both panels of Fig. 11.51 the curve of this orbit lies on the 
outside. Hence the resonant islands first appeared just inside this curve. As 
the potential flattened more, fi^/fiz dropped significantly below unity for the 
planar orbit, so the resonance condition was satisfied by orbits with non-zero 
and the islands moved inwards. Orbits whose curves lay in the path of a moving 
island did one of two things: (a) they were trapped into the island, or (b) they 
abruptly increased their radial actions so that their curves went round the far 
side of the island. Which of these two outcomes happened in an individual case 
depended on the precise orbital phase of the star when the potential achieved a 
particular flattening, but it is most useful to average over phases and to consider 
the outcomes to occur with probabilities Pa or P^. The magnitudes of Pa and 
Pb depend of the relative speed with which the island increased its area and 
moved: if it simply grew. Pa = 1, and if it moved without growing Pb = 1. 
These results follow from Liouville's theorem that phase-space density cannot 
increase. 

Let's imagine that after a period of stationary growth, the island moved 
inwards without growing, and then became stationary while it shrank. In this 
case it would have swept up stars with large and small Jz and released these 
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stars into orbits with smaller Jr and lar ger J^. In other words, it will have 
turned radial motion into vertical motion. Sridhar fc Touma (Il996l) have called 



this process "levitation" . Conversely, the moving island will have reduced the 
vertical motions and increased the radial motions of any stars it found in its 
path through action space. Thus resonances stirr the contents of phase space. 
Levitation is a lovely idea but it's not clear that it is of practical importance. 
We shall see below that in a disc analogous scattering by resonances is very 
important. 

1.2.4 From order to chaos 

If we set z = and express the triaxial potential (|1.39|) in cylindrical coordinates, 
it becomes 

0) = ^vl In [rg + ii?2(^-2 ^ ^) _ ii?2(g-2 _ 1) cos 2^] . (1.40) 
Consider now the related potential (jBinnevl . Il982l) 



ii?2((7-2 + 1) _ (^ii?2(q-2 - 1) + — j cos2</) 



(1.41) 



Fig. ll.lOl shows surfaces of section for motion in this potential when rg — 0.1, 
q — 0.9 and Rc — 6 (left) and i?o = 4 (right). Comparing the left panel of this 
figure with the left panel of Fig. II. 61 we see that the extra term in the logarithm 
has greatly increased the number of resonant islands, and in the right panel 
we see that it can introduce a qualitatively new feature: many points seem 
now to be scattered at random rather than confined to curves. Confinement to 
curves is an indication that the orbits admit an isolating integral in addition 
to energy, and are in consequence quasiperiodic. Conversely, when the points 
are not confined to curves, the orbits lack an additional integral and are not 
quasiperiodic. We say these orbits are irregular or chaotic. 

The increase in the number of resonant islands is probably the cause of the 
emergence of chaos. In the surface of section, islands form chains, with the 
curves of some non-trapped orbits running on one side of the chain, and those 
of others running on the other side of the chain - Fig. ll.4l mav help the reader to 
picture this situation. The curves of the most nearly trapped orbits on each side 
come very close where the islands touch. The tiniest perturbation can cause a 
star on one of these orbits to swap from one side of the chain to the other. If a 
star makes such changes, its orbit ceases to be quasiperiodic. 

If there are several chains of islands in the surface of section, and the islands 
of two chains almost touch, a star can make two such swaps, moving from, say, 
inside chain 1 to outside chain 3. In this way a star can move stochastically 
through a significant region of phase space. This is probably what happens in 

Fig.oni 
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Figure 1.10: Surfaces of section for the potential ()1.4ip with q = 0.8 and ro = 
0.1. The left panel is for — 6 and the right panel is for Rc — 4. 




Figure 1.11: A fast encounter between two stars at speed v and impact param- 
eter b. The force perpendicular to the relative velocity is ~ Gm,im2/b^ and it 
acts for a time ^ 2b /t. 



1.3 Fluctuations 

In previous sections we investigated the orbits of stars in a smooth, time- 
independent model of the galaxy's gravitational potential. In reality the po- 
tential contains time-dependent features and in this section we investigate how 
these features drive evolution. 

A fundamental result is obtained by multiplying the equation of motion 
p = — V$ by p and rearranging the result to 

Thus stars change their energies if and only if the potential is time-dependent. 
Fluctuations in the potential enable stars to exchange energy. 



1.3.1 Two-body scattering 

The most obvious source of fluctuations is the moving gravitational potentials 
of individual stars. When stars of mass mi and m2 pass each other at speed 
V and impact parameter b (Fig. Il.ll|) the effect is an exchange of momentum 
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along the line that is perpendicular to the mutual velocity and has magnitude 
~ 2Gmim2/bv. So the encounter adds to the velocity of mi a velocity i5vi of 
magnitude ~ 2Gm2/bv. The direction of these increments is random, so we add 
them in quadrature. The rate of such encounters is ~ 2Tmvbdb, where n is the 
number density of stars, so the rate of change of I'^'^iP is 

-^|JviP = 2.nyd66(^^j ^y-. (1.43) 

The integral diverges at both ends of the range of integration. The divergence 
at small b is an artifact that can be traced to our use of 2Gm2/bv as the magni- 
tude of the velocity change in an enc ounter: an accurate calcul ation shows that 
the velocity change never exceeds v ( Binnev Sz Tremaind . 20081 eq. 3.53a). The 



divergence at large b is real, and indicates that encounters with impact param- 
eters that are on the order of the size of the system dominate. Physically, what 
this means is that the dominant source of fluctuations is Poisson fluctuations 
in the number of stars in substantial parts of the system. The mass inside a 
volume of radius r will fluctuate by SM ^ M/\/N, where N = ^-Kr^n is the 
number of stars in this volume. Just outside this volume the gravitational field 
will fluctuate by 



G5M GM Gm2 „ ,4 



ViV = Gm2W|^- (1.44) 



This fluctuation acts for a time ^ r/u so it changes the velocity of any star by 

Gto2 



|^v| -^^^irnr, (1.45) 

which grows with r, thus confirming that large-scale fiuctuations are the most 
effective. If we accept that the dominant fluctuations are those involving half 
the system, so r is about half the system size R, we conclude that the number 
of half-crossing times r/v required to change v by of order itself is 



1 



|<5v|2 (Gm2)2 ImiR 
The virial theorem implies that ~ GM/ R, so 

t2B f GM Y 1 (M/ma) 



Across \RGm2 J ^TmR ^nnR^ 



(1.46) 



N. (1.47) 



For a galaxy this number of half-crossing times is many times the age of the 
Universe, so the fiuctuations associated with the motions of individual stars are 
unimportant. But for a globular cluster, which has N < 10^, i? ^ 3pc and 
V ^ 6 km see" so r/v ^ 0.25 Myr, v changes by order itself in ^ 6 Gyr so the 
process is significant. In an open cluster the process is even more important. 
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Two-body interactions randomise the distribution of stars in phase space and 
thus drive the system towards thermal equilibrium. No such e quilibrium is possi- 
ble fo r a stellar system that is only confined by its own gravity (jBinnev fc Tremaind . 
20081 §4.10). But we can understand the impact of two-body interactions by 
considering the consequences of trying to reach thermal equilibrium. 

The speed required to escape from a stellar system is never much larger 
than the system's characteristic velocity dispersion a - one can easily show 
from the virial theorem that the mass-weighted rms of the local escape speed 
is only twice the mass-weighted rms velocity dispersion: (wo)^^^ = 2(cr^)^/^. 
Consequently, the velocity distribution is always distinctly non-Gaussian. Two- 
body scattering drives the velocity distribution towards Gaussianity, so it is 
constantly trying to repopulate the missing tail of the velocity distribution at 
V ^ 2a. Stars scattered into this domain are free and leave the system, to the 
system loses mass by evaporation on the two-body timescale. 

In thermal equilibrium there would be equipartition between the particles. 
So massive stars would have a smaller velocity dispersion than low-mass stars. 
Consequently, two-body interactions are constantly transferring energy from 
more massive to less massive stars, with the consequence that the massive stars 
sink towards the centre of the system: two-body scattering drives mass segre- 
gation. 

In thermal equilibrium all parts of a body have the same temperature. In 
a self-gravitating system there is a tendency for the centre to be hotter than 
the outside, if only because the escape speed decreases outwards. So two-body 
interactions tend to transfer energy outwards from the core to the envelope. 
By the virial theorem, a self-gravitating system that loses energy contracts and 
gets hotter, while one that gains energy expands and becomes cooler. So the 
conduction of heat from the core to the envelope increases the difference in 
temperature between the two parts of the system and accelerates the heat flow. 
The upshot is the gravithermal catastrophe in which the core contracts in 
both size and mass until it contains only a few stars. 

The point to note about evaporation, mass-segregation and the gravithermal 
catastrophe is they are all consequences of fluctuations in the gravitational field 
driving the system towards an unattainable thermal equilibrium. In star clus- 
ters fluctuations associated with individual stars are sufHcient to generate these 
effects on astronomically interesting timescales. In galaxies they are not, but 
evaporation and the gravithermal catastrophe will be driven by whatever fluc- 
tuations do occur, while equipartition won't be because it depends on stars of 
different masses experiencing different fluctuations. Sources of significant fluctu- 
ations in galaxies include giant molecular clouds, spiral arms, satellite galaxies, 
and high-speed encounters with other comparable galaxies. 

1.3.2 Orbit-averaged Fokker-Planck equation 

In this section we develop a general framework for handling the impact of fluc- 
tuations. The general idea is that, by the strong Jeans theorem f §1.1.2p the 
galaxy's distribution function is at all times a function f{3,t) of the actions. 
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Fluctuations (and resonances) cause this function to evolve by causing innumer- 
able small changes i5J in the actions of individual stars. Let P(J, A)d^ A St be 
the probability that in time 5t a star with actions J is scattered to the action- 
space volume d^ A centred on J + A. The number of stars in the action-space 
volume d'^J is (27r)'^/(J, i)d'^J, so the number of stars leaving this volume in St 
is 

{27rf f {3, t)d^ 3 St J d^AP{3,A). (1.48) 
Similarly, the number of stars that are scattered into this volume is 

{2TT)^d^3St J d^A /(J - A, t)P{3 - A, A). (1.49) 

Hence the rate of change of the distribution function is 

J d'A [/(J - A, i)P(J - A, A) - /(J, t)P{3, A)] . (1.50) 

Since scattering events change actions only slightly, P{3, A) is appreciable only 
for |A| ^ |J|. So we can truncate after just a few terms the Taylor series 
expansion in J of the product /(J, t)P{3, A): 

/(J - A, t)P{3 - A, A) = /(J, t)P{3, A) 



dJ,, + 2 QJ^QJ^ + (1-51) 



Substituting the first three terms on the right side of this expression into equa- 
tion (|1.50p and cancelling terms, we obtain 



A,(J) =y"d3AA,P(J,A) and A2^.(J) = ^ d^A A,A,F(J, A). (1.53) 

Equation (|1.52l) is the orbit-averaged Fokker-Planck equation. It states 
that the rate of change of the distribution function is minus the divergence of 
the flux F of stars in action space, and we have an expression for that flux in 
terms of the diffusion coefficients defined by equations (|1.53p . The latter are 
simply the expectation value and the variance of the probability distribution of 
changes in actions per unit time. 

The diffusion coefficients reflect the physics of whatever is responsible for 
causing the fluctuations. In some circumstances, for example in a star cluster, 
the fluctuations will be approximately thermal in nature, with temperature T. 
Then the principle of detailed balance requires that the stellar flux vanish when 
the objects being scattered are in thermal equilibrium with the fluctuations. 
That is, F = for 

/(J) = constant x e"^/''^, (1.54) 
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where H{3) is the Hamihonian. In this case we have 



so F = implies that 



= / 



(1.56) 



Clearly the square bracket must vanish, so we obtain an express ion for the first- 
order diffusion coefficient in terms of the second-order coefficient ([Binney fc Lacev 
Il988l) 

"-i'^-I^A?. . (1.57) 



2 



dJj kT 



This expression is useful because it enables us to obtain the first-ord er d iffu- 
sion coefficients from the second-order diffusion coefficients Afj, and, 

while Afj can be obtained from first-order perturbation theory (see below), a 
direct calculation of A^ requires second-order perturbation theory. 

The diffusion coefficients are conveniently calculated by expanding the po- 
tential in angle-action coordinates 

$(x,t) = $o(x) + $i(x,t) ^ 'i'o + ^^^{3, t) cos{n ■ e + ijjr,). (1.58) 

n 

where $o(x) is the potential of the underlying Hamiltonian Ho{J) and is the 
fluctuating part of the potential. Hamilton's equation of motion for J is 

j - = = E ^) • ^ + ^n)- (1-59) 

n 

To get a random change in J, we need to integrate this equation of motion for 
a time T that is longer than the auto-correlation time of the fluctuations. We 
do this by expanding the variables in powers of $i/$o: 

J(t) = Jo + Ai(<) + A2(i) + --- and e{t) = Oq + Qot + Oiit) + ■ ■ ■ (1.60) 

Ai is obtained by integrating equation (|1.59p along the unperturbed orbit 

Ai(r) = Vn /" dt$n(J,i)sin(n-6/ + Vn) 
„ -Jo^ 

= Vn /" di $„( J, i)sin[n-(6>o + m) (1.61) 
„ -^0 

To obtain the second-order diffusion coefficient we multiply this equation by 
itself and average over initial phases Oq. After reordering the integrals so 9o is 
integrated over first, we find that the innermost integral is 

(27r)"3 j d^0^ sin[n • {Oo + flot) + v^„] sin[n' • (6»o + flot') + Vn']- (1-62) 
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Using 2 sin ^ sin B = cos(A — B) ~ cos{A + B) and that the integral of any 
cosine that depends on will vanish, we conclude that the innermost integral 
vanishes unless n' = n0 when it's equal to ^ cos[n • ^^{t — <')]. Hence 

(AhAi, (T)) = iVn,nj / dt / di' $„(J, t)$„(J, t') cos[n • f2o(t - t')]- 
„ Jo Jo 

(1.63) 

Next we take the ensemble average over the fluctuations that are represented by 
$n- We assume that they are a stationary random process so the autocorrelation 
of $n(J, depends only on the time lag t — t': 



$n(J, i)$„(J, t') = c„(J, t-t'). (1.64) 
with this assumption we have 

Ah Aij (T) = iVn.rij / dt / dt' c„(J, t - t') cos[n • r2o(t - t')] 
„ Jo Jo 

dt;c„(J,w) cos(n • Jlo^') / du (1.65) 

'T J\v\ 

I duc„(J,w) cos(n • r2ow)(r — jwl), 

J-T 



where in the second line we have introduced new coordinates u = t + t' and 
V = t — t' . Given that we want T to be bigger than the autocorrelation time 
of the fluctuations, we have that whenever Cn{3,v) is non- negligible, \v\ <C T, 
so term in the integrand that's proportional to \v\ can be neglected, leaving 
a result that's proportional to T. The diffusion coefficient is the coefficient of 
proportionality, so 

A|^= i^n,nj2-„(J,n-Oo), (1.66) 

n 

where Cn(J,^) is the power spectrum of the fluctuations: 

Cn(J,w)= / dw Cn(J, w) COS(CJW) = / du <I>n(J, t)$n(J, t — w) COs(aJv). 
J-T J-T 

(1.67) 

The bottom line of this result is that the ability of a star to diffuse through 
phase space hinges on whether the fluctuations contain power at one of the star's 
natural frequencies n • rjp. In particular, if the fluctuations are periodic in time, 
for example because they arise from a normal mode of the system, they will drive 
diffusion only of stars that resonate with them. In practice periodic fluctuations 
will simply depopulate narrow regions of phase space: stars for which n • rig is 
equal to the frequency of the fluctuation will be scattered to new actions and 



^ Since we are using cosine series, we need sum over only half of n space so n' will never 
equal — n. 
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then ce ase to be resonant becaus e fundamental frequencies are functions of the 



actions. ISehwood fc KahnI (|l99lh find evidence for such action-space "grooves" 
in numerical simulations of stellar discs and show that they can generate new 
spiral features, which in their turn generate other grooves. 

1.3.3 Heating of the solar neighbourhood 

Fig. 11.121 shows the radial, vertical and azimuthal velocity dispersions of groups 
of nearby stars with accurate space velocities as a function of the colour of the 
stars. Blue stars are plotted on the left and red stars on the right, so all three 
components of velocity dispersion increase from blue to red. Blue stars are 
massive and short-lived, so in the blue bins all stars are quite young, while red 
stars live longer than the age of the galaxy, so in the red bins we have stars of 
all ages, but with a bias to old stars because the star-formation rate was higher 
in the past than it is now. So the variation of velocity dispersion with colour 
indicates that the random velocities of stars increase over time. From these data 
and isochrones one can deduce how velocity dispersion i ncreases with age, and 
the conclusion is that a ~ (jAumer fc Binnevl . l2009| ). 



It's instructive to infer from this result how the diffusion coefficients must 
scale with |J|. We make two simplifying assumptions: (i) that the dominant 
scatterers are much more massive than stars, and (ii) that the velocity disper- 
sions of groups of stars scale with the mean actions in the group as 



Cr oc \/(Jr) and <7z oc ^yjjz). (1.68) 

These relations are exact in the epicycle approximation, in which the radial and 
vertical oscillations of stars are harmonic, so for example Jr = Eji/n^ Since 

^ Quite generally we have that QrJr is equal to the time-averaged value of along any 
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scattering must be dominated by giant molecular clouds and spiral arms, the 
assumption of massive scatterers will be a good one. In thermal equilibrium 
with such massive bodies, stars would have velocity dispersions that are larger 
than those of the clouds and arms 7 km see" ^) by the square root of the 
ratio of masses, so the stars' velocity dispersion would be > lOOOkmsec"^. 
Consequently, we can use equation (ll.57|) in the limit of infinite temperature|f| 
when the Fokker-Planck equation simplifies to 

Stars are born on orbits that have non-negligible angular momenta = 
but small values of Jr and J^. Consequently, a young population is initially 
distributed in action space along the axis, and diffusion of this population is 
predominantly away from this line, towards larger values of Jr and . For this 
reason we neglect derivatives with respect to in equation (|1.69p . 

In problems involving the ordinary diffusion equation, a key solution is the 
Green's function exp(— x^/2<)/(27r<)^/^, which describes the spatial distribution 
at time t of particles injected at a; = at time t = 0. Analogously, we seek a 
Green's function of the form 

/ = i-2a/„(X) where X = ^. (1.70) 

In this solution the mean value of | J| will increase with time as t"-, and the power 
of t multiplying /o ensures that the total number of stars / dLz J dJrdJz f is 
conserved as stars diffuse from the axis. Suppose Af^ scales such that Afj{k3) = 

A:''A|'(J). Then putting k = t-'' we have A|"(X) = i-"''A|^(J). Evaluating 
both sides of equation (|1.69p with these assumptions yields 

dfo\ i,ab-ia 9 f-rT/-v\'^fo 



^2«/o + .X.^J — ^A^.(X)^J. (1.71) 

This equation can be valid at all times only if 2a + 1 = 4a — afe, so = 2 — 1/a. 
Consequently, the empirical result ( J^) ^ cr^ ^ t^/^ implies a c± | a nd 6—5- 
The scaling (7^ ^ t^^"^, which has been advocated by IWielenI (Il977h and 



several subsequent authors, implies a = & = 1. A simple argument shows that 
it is implausible for the diffusion coefficients to grow so rapidly with |J|. In 
the epicycle approximation, J^ differs from the epicycle energy E]^ only by 
the (constant) epicycle frequency, so ^ AEn. = v • (5v, where (5v is the 
projection into the equatorial plane of the change in a star's velocity as a result 
of a scattering event. Hence (A^) ^ |J| implies 

Eji ^ {(AEnf) ^ ((v • 5v)2) ^ (iJflJJvp). (1.72) 



orbit. 

^See Appendix B of lBinnev &: LacevI | |198S|) for a rigorous justification of this step. 
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That is, CTr ^ t^^^ implies that |i5v| is independent of |v|. However, gravitational 
scattering always causes the momentum change 5v to decrease with increasing 
speed because the gravitational force is independent of speed and the time for 
which it acts decreases as 1/1 v|. 



Can we derive Afj{k3) ~ fc^/^Afj(J) from physics? iBinnev fc LacevI (|l988l ) 
show that this sca ling is predicted by the model of cloud-star scattering that 
was introduced by [Spit zer fc Scwarz schildl (|l953l) . However, this model is de- 
fective in two respects: (i) it assumes that the relative velocity with which a 
star encounters a cloud is dominated by epicycle motion rather than differen- 
tial rotation, and, more seriously, (ii) it assumes that stars are confined to the 
equatorial plane. In reality as a star ages it oscillates with increasing amplitude 
and period perpendicular to the plane, and these oscillations decrease its prob- 
ability of being scattered by a cloud. Consequently, when this effect is taken 
int o account, Afj(J) increases with |J| more slowly than as |J|^^^. 

Binnev fc LacevI (fToSS) show that three-dimensional scattering by molecular 



clouds generates a tensor of diffusion coefficients Afj which is highly anisotropic. 
The consequence of this anisotropy is that we exp ect az/o ' r -8, which is 
significantly larger than the observed value, ^ 0.6. S ellwoodI (|2008r ) argues that 
the discrepancy arises from the erroneous assumption of an isotropic distribution 
of encounters: as in two-body scattering, distant encounters are important, and 
since both stars and clouds lie within the disc, distant encounters are dominated 
by the velocity components that lie within the plane and do not change Jz- 

Thus it seems that scattering of stars by giant molecular clouds may set 
the ratio of the vertical and horizontal velocity dispersions of disc stars. While 
star-cloud scattering makes a significant contribution to the secular increase 
in the velocity dispersions of stars, it probably cannot account fully for the 
data because its effectiveness declines rapidly with increasing velocity dispersion 
and thus cannot account for the numbers of stars with radial dispersions > 
30 km sec" ^. 



1.4 Spiral structure 

Thin galactic discs are extremely prone to generating spiral structure. Many 
manifestations of spirality arise from gas-dynamics rather than stellar dynamics 
- for example the chains of blue stars so evident in blue and UV exposures, 
and the spiral distributions of HI and CO detected at radio frequencies. In 
such cases the old stellar disc is believed to carry a large-amplitude spiral also, 
but it is not easy to detect these a stellar spirals. It is most evident in near 
IR photometry, which is dominated by a combination of stars near the to p 



of the giant branch and low-mass main-sequence stars (|Rix &: Zariskvl 119951 ). 
The former are a measure of recent star formation so they tell us more about 
gas than stellar dynamics, but the latter conta in most of t he m ass of the disc. 



so are a window into the key dynamics. iRix fc Zariskvl (1995) find that the 



2.2 /im surface brightnesses of spiral galaxies carry spiral structures that have 
amplitudes of order unity, and that the fluctuation in the surface density of stars 
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Figure 1.13: Energy versus angular momentum for planar orbits in an axisym- 
metric potential - a "Lindblad diagram". No orbits lie in the shaded area, 
which is bounded by the points of circular orbits. A potential that is stationary 
in a rotating frame moves stars along lines with slope dE/dLz = flp- (From 
Sellwood & Binney 2002) 

is probably nearly as large. Radio-frequency spectral lines and the Ha line show 
that spiral disturbances are associated with streaming velocities ~ 7 km see" ^. 

Spiral structure proves to be an intrinsically non-linear phenomenon, and as 
a consequence our understanding of it is still frustratingly incomplete. It's raison 
d'etre is, however, clear: it is the principal means by which galaxies transport 
angular momentum outwards, which enables them to increase their entropy - 
i.e., heat their discs. We start our study of spiral structure by examining how 
this heating comes about. 

1.4.1 Secular evolution driven by spiral structure 

Let's assume that spiral structure is a nearly stationary pattern that rotates at 
some fixed angular speed ftp. In this case, when we write the potential of the 
spiral structure in angle-action variables (eq. I1.58|) . the expansion coefficients 
$n(J,i) will contain only multiples of fip in their temporal Fourier transforms. 
Hence the power spectrum of the potential Cn(J,i^), which appears in equation 
(jl.66l) for the diffusion coefficients, will be non-zero only when w is equal to one 
of these frequencies, so the spiral will enable only a minority of stars to diffuse. 
We now examine the impact of the spiral on the minority stars that resonate 
with it. 

If we work in the frame of reference that rotates at frequency fip, the motion 
of each star is governed by a time-independent Hamiltonian, the numerical value 
of which, the Jacobi constant, is an isolating integral. In terms of the energy 
E of motion in the non-rotating frame, the rotating-frame Hamiltonian is 



H — E — ^IpLz- 



(1.73) 



James Binney Dynamics of secular evolution 



29 



Since H is an integral, dH = and changes in E and caused by the spiral 
satisfy 

dE^QpdL^. (1.74) 

Fig. 11.131 is a plot of E versus Lz , and equation (|1.74p states that in this figure 
a steadily rotating spiral moves stars on lines of slope Jlp. The physically ac- 
cessible region is bounded below by the locus of circular orbits, which are the 
orbits with the largest value of for each given E, so there are no orbits in the 
shaded region below this boundary. The slope of the boundary, {dE/dLz)j^=o, 
is the circular frequency Vl{Lz)- Clearly at the corotation resonance (CR), 
where ri(L^) = fip, the spiral scatters stars from one circular orbit to another. 
Elsewhere, the spiral scatters stars away from the boundary, to places where the 
energy exceeds that of the circular orbit of the given value of L^, and the ad- 
ditional energy will be invested in epicyclic motion. Inside the CR, the angular 
momenta of stars must be reduced, while outside the CR it must be increased. 
Thus the spiral must move angular momentum outwards. 

We have seen that significant shifts in actions only occur at resonances, 
where n • fJo — mSlp, where m is the number of arms that the spiral has be- 
cause the time dependence of the potential is cx cos{m^pt + ip). Besides the 
CR [n — (0, 771, 0)], the two most important resonances are the inner Lindblad 
resonance (ILR), where n = ( l,m,0), and the outer Lindblad resonance 
(OLR), where n = (1, to, 0). At the ILR the Doppler-shifted frequency at which 
a star perceives the spiral is m(Jl — Vl-p) and this coincides with its radial fre- 
quency, rir, while at the OLR the perceived frequency of the spiral is m(yi-p — 17) 
and this again coincides with 17^. We have shown that the spiral absorbs at 
the ILR and emits it at the OLR. At both places changes in heat the disc. 

At the CR the change in angular momentum can have either sign, but the 
star simply moves from one circular orbit to another, so the disc is not heated. 
In fact these shifts of stars are so inconspicuous that for decades they were 
overlooked. They are astronomically important, however, because in galactic 
discs metal abundances generally decrease outwards, so the radial migration 
of stars at the CR can be detected if metallicities are measured. Specifically, 
radial migration ensures that at each radius there are stars of the same age but 
dif fering metal abunda n ces be cause they formed at different radii. 



Sellwood fc BinnevI (|2002| ) explored these phenomena with N-body simula- 



tions of discs. In one experiment they introduced an action-space groove into 
the initial conditions to generate an isolated, transient spiral feature. The left 
panel of Fig. 11.141 shows the distribution of ensuing changes in versus initial 
Lz- Vertical lines mark the locations of the CR and Lindblad resonances for 
the measured pattern speed. Stars interior to the CR gain and transfer to 
outside the CR, while those outside the CR lose and move inwards. Thus 
these stars swap places. The right panel explains how this is done by plotting 
orbits in the ((/), R) plane when a steady two-armed perturbation is imposed. 
There are islands formed by orbits that are trapped at the resonance, and wavy 
lines of orbits that continue to circulate. Orbits forming the island constantly 
move from inside to outside the CR and back again. When a spiral potential 
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Figure 1.14: Left: the distribution of changes in angular momentum amongst 
the first 20% of stars in an N-body simulation when they are ordered by initial 
epicycle energy. The vertical lines show the angular momenta corresponding 
to the ILR, CR and OLR of the spiral pattern within the simulation. The full 
curve shows the mean of the distribution and the shaded region is bounded 
by its 20th and 80th percentiles. The dashed line has a slope of —2. The right 
panel shows the response of initially circular orbits to a transient spiral potential 
(after Sellwood & Binney 2002). 
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Figure 1.15: Each panel shows the distribution of the final guiding-centre radii 
of the stars in a disc simulation whose initial guiding-centre radii lay within the 
region between the dotted vertical lines. The disc had a flat rotation curve, 
Q = 1.5, and half of the radial force was provided by a fixed halo. The duration 
of the simulation was ~ 4 Gyr. (From Sellwood & Binney 2002). 



emerges, it creates islands which grow with the potential by sweeping up orbits 
from the wavy regions. When the potential fades, the islands shrink and the 
trapped orbits are released on each side. A star that was on an inner wavy 
orbit at entrapment may move on a trapped orbit from inside the CR to outside 
the CR before being released to a wavy orbit outside the CR. The maximum 
distance stars can move their guiding centres is set by the largest extent of the 
islands. Stars that start far inside the CR are released far outside the CR. The 
tendency for the populated regions in the left panel of Fig. ll.l4 | to slope from top 
left to bottom right with gradient —2 confirms this picture. ISellwood fc Binnev 
called this process of swapping places around CR churning. 
In another experiment Sellwood & Binney allowed a disc to evolve for ~ 
5 Gyr without seeding any spiral structure. Irregular and quite weak spiral 
structure emerged as the disc was evolved. Fig. 11.151 shows histograms of the 
final radii of stars that all started from the narrow radial bands marked. Stars 
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that started at the current radius of the Sun finished at radii that are often 1.2 — 
2kpc from Rq. From the strength of spiral structure seen in NIR photometry, 
Sellwood & Binney estimated that over the Hubble time stars will typically 
migrate ~ 2 kpc from their birth radii. 

1.4.2 Propagation of spiral waves 

If we disturb the surface of a pond with a stone, the water molecules hit by the 
stone disturb water molecules slightly further away, which in their turn disturb 
their neighbours, and in this way much of the energy of the original impact 
is carried over the surface of the pond to be dissipated as the waves hit the 
pond's shore. If we disturb the stars in some region of a disc by perturbing 
the gravitational potential in their neighbourhood, the orbits of these stars will 
change, and a short time later their contribution to the disc's density distribu- 
tion will change, which will modify the disc's potential. This modified potential 
will disturb the orbits of other stars, and in this way a way of disturbance will 
propagate over the disc. 

A significant difference between the pond and the disc is that in the pond 
water molecules excite their neighbours through prcssiire: they simply push 
on molecules they touch; the interaction is local. In the disc stars disturb other 
stars by modifying the gravitational field, which propagates at the speed of light. 
The latter is effectively infinite: we are dealing with action at a distance, so the 
physics is inherently non-local. If we take this non-locality seriously, we have 
to compute globally and proceed straight to the determination of the normal 
modes of the entire disc. 

There is, however, a special case in which the gravitational interaction is 
effectively local. This is the case of tightly wound waves: as one moves away 
from a density wave in the disc, the latter's gravitational field decays on the 
lengthscale of a wavelength because, by virtue of the long-range nature of grav- 
ity, the positive and negative contributions to the field at the point of observation 
from peaks and troughs soon cancel rather precisely. If waves in a disc have a 
wavelength that is small compared to the local radius, their gravitational field 
is highly localised. For such waves we can derive a dispersion relation and thus 
obtain the group velocity, etc. 

We now set up the equations whose solution yields the normal modes of a 
stellar disc. Then instead of solving them we use the tight- winding approxi- 
mation to derive the dispersion relation of tightly wound spiral waves. Finally 
we use this dispersion relation to gain insight into the global dynamics of self- 
gravitating discs. 

We start by finding how the DF /(J) is changed by a perturbing potential 
5<i>(x). The governing equation is the collisionless Boltzmann equation 




(1.75) 



where [..] denotes the Poisson bracket. Writing f{J,d,t) = /o(J) + fi{J,0,t) 
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and H{J,6,t) ~ Ho{J) + Hi{3, 9, t), we have to first order 

fi + rf, fi_ dh dH, dfo 

— - [Ho,h] + [H,Jo] - . — + — .—. 



(1.76) 



We note that OHq /dO = ^Iq and Fourier expand Hi and /i 
/i(J,0,t)=.^5/„(J)e'("-«-*) 

n 

Til (J, e,t)=Y, 5<i>„(J)e'("-«-'^*), (1.77) 

n 

where by virtue of the time-translation invariance of equation (|1.76l) we have 
assumed that all perturbed quantities have the same frequency, lu. Using these 
expansions in equation (|1.76l) we equate coefficients of e'("'^~'^*) on each side to 
obtain 

SU - (1.78) 
n • S2 — 

The frequencies w of the system's normal modes are determined by the 
requirement that $i is related by Poisson's equation to the density fluctuation 
implied by /i. This requirement reads 

^ ^ (^<i>„e--) = e'- / d^P A = E / 

n •' n •' 

(1-79) 

This equation is homogeneous in so we expect it to have non-trivial solutions 
only for particular values of w, the frequencies of normal modes. Finding the 
normal-mode frequencies is hard because the equation involves in an essential 
way two systems of phase-space coordinates: (x, p) and (0, J). 

Although equatio n (|1.79p can be tackled (see for example iKalnais . 1977 : 
Read fc Evamil998l ). we will pursue a simpler course. We restrict ourselves to 



razor-thin discs, which have a four-dimensional phase space, and use the tight- 
winding approximation. We assume that the disturbance is a tightly wound 
spiral wave, so 

$i(i?,0,t) = ee'('=-«+'"'^-"*) (1.80) 

with kR ^ 1 for trailing waves and kR <C — 1 for leading waves. It is straight- 
forward to show from Poisson's equation that the corresponding surface density 
is 

Si--;^*!. (1.81) 

We adopt the epicycle approximation (i jl.l.Sp . in which the real-space coordi- 
nates are related to angle-action coordinates by 

R^R^ + acos9r, = 6*0 -t- ^ sin 6*^ , (1.82) 

Rg 

where 



2 7 

— and 7 = 2n/K. (1.83) 

K 
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Then the Fourier decomposition of $i is 



:ee"=^«exp(ifcacos6l^)e'"''*exp i— ^sin^^ e""^*. (1.84) 



We now combine the two exponentials of circular functi ons into a single expo- 
nentia l of a single cosine and use equation (8.511.4) of iGradshtevn fc Rvzhik 



;o express this as a sum over Bessel functions J7/. Specifically 



exp 



1 I ka cos Or H — sm ( 



— exp [ia/C sin((?r + a)] 

= ^ J/(/Ca)e"^'''-+"^ (1.85) 



where a is the pitch angle and /C is the total wavenumber of the spirals: 



/ 7717 \ / TJlP"^'^ 

a{J^) EE arctan {^-^ J and /C( J^) = J + . (1.86) 

Using equation (|1.85p to rewrite equation (|1.84p we obtain 

CO 

$l(6>,J,t) = e e'^kR,+lc.) j^(^j^^y{W^+m9^-ut) ^ (^g^^ 

Equation (|1.87p tells us what (5$n(J) is for n = (/, m, 0): 

'^'i>(;,m.o) =ee'('=^-+'")ji(/Ca). (1.88) 
Using this in equation (|1.78p we obtain the change in the DF caused by $i: 

SU = "■^/o/^J ,eH^H,+Mj-,(^^) [„ = (^,^,0)]. (1.89) 
n • S2 — 

The change in the surface density is Ei = / d^p/i. Since our expression for 
/i uses angle-action coordinates rather than (x, p) coordinates, we use a trick 
based on the fact that d^a;d^p — d^^d^J because both systems are canonical: 

Si(i?>,i)= J d^pfi^^J d^P J AR'R' j d(b'6{4>~cl>')6{R~R')h 

= ^J d^Jd^eSicf) - (l)')S{R - R')fi (1.90) 

f cfJd^e S{(l> ^ ^')SiR ~ R')J2^fne'^"'^-^'l 

Using equation (|1.82p for and i?', and equation (11.89^ for ^/„, this becomes 

Y.i{R,(f),t) = ^j d^ Jd^e (5 {j^-^'p-'Y S{R-Rg-a cos Or) 
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n • r2 I 



JiiICa) 



(1.91) 



The two Dirac delta-functions enable us to carry out the integrals over 0^ and 
J^. This done every occurrence of Rg{J^) (including those in the definitions 
of K, a, etc.) should, strictly, be replaced by i? — acos6'r.. However, the tight- 
winding approximation allows us to neglect the small difference between i?g and 
R except when it occurs in the argument of an exponential multiplied by the 
large number k. With the aid of this approximation we obtain 

l]i(i?,0,i)^|e'('=«+"'^-*)^ V e''" / dJ.J-KM "'y^ 



l=-c 



X / d9r exp 



i (j,6r — sinOr ~ ka cos 9r^~^ [n = (^, m, 0)]. (1-92) 



It is simple to show that dJ0/di?g = dL^/di?g = R^k/j. For /o we adopt 



f (1) - '^^0 -^,J,/a^ 



(1.93) 



which with the epicycle approximation (eq. I1.82p yiel ds the Schwarzschild veloc- 
ity dist r ibutio n with radial dispersion a (e.g. ^54.4.3 [Binnev fc Tremaine . 2008 : 
Binnevl . [2010l ). Finally, we use equation (|1.85p to express the exponential of si- 



nusoids in the last line of equation (|1.92l) as a sum over Bessel functions. Then 
we can evaluate the integral over Or to obtain 



El(i?,(^,i): 



0, gi(fci?+m0-Ljt) 



_ ^^^0 i(feH+m0-wt) 



E 

— oo 
oo 

E 



-I 



Ik + mfl — u) 



Ik + mfl — u! 



dJr\Ji{ICa)\\-^J^/^" 
(1.94) 



where equation (6.615) of iGradshtevn fc RvzhikI (|1965I ) has been used to evalu- 
ate the integral over J^, Ii is a modified Bessel function, and 



(1.95) 



Two more definitions and the identity Ii{z) — I-i{z) enable us to write equation 
p.94p in the neater form 



Ei{R,(f>,t)- 



k2(1-s2) 



where 



uj — mfl 



Y 1 



hix) 



2 /;2- 



1=1 



S^/l 



(1- 



(1.97) 
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The final step is to require that the value of Si from equation (|1.96p agrees 
with that given by equation (|1.81l) . Eliminating Ei between these equations and 
approximating /C by |fc| (see ea. ll.86p . we obtain the Lin-Shu-Kalnajs dispersion 
relation for tightly wound spiral waves: 

|fc| 2\ ("^ ^ jnQ) 



•^(■s,x) = (1 - s ) = 1 — , where fccrit = 2^^^^ ' 



For an m-armed spiral with pattern speed Vlp, uj = milp so 



(1.99) 



which rises from —1 at the ILR through zero at the CR to 1 at the OLR. Hence 
from equation (|1.98p kT vanishes at the Lindblad resonances. One finds that 
T behaves roughly as k(\ — bk) with b > , so the left side of the dispersion 
relation peaks for some fc, and Toomrd (1964) showed that if the disc is stable to 



axisymmetric disturbances this peak value is smaller than unity, so solutions for 
k cannot be found for a range of small values of s^. Thus waves are forbidden 
in a zone around the CR as well as inside the ILR and outside the CR. In the 
permitted zones two values of k can be found for given s, the values approaching 
one another as s approaches the forbidden zone around the CR. Thus there are 
two branches to the dispersion relation, and the branches merge at the edge of 
the CR z one. 

When lToom^ (|l969l) determined the group velocity of waves from the dis- 



persion relation, he found that short-leading waves propagate outwards from the 
ILR. At the edge of the forbidden region around the CR these waves transfer 
to the long-leading branch and propagate back towards the ILR. As the waves 
approach the ILR, k decreases and the validity of the tight-winding approxima- 
tion becomes questionable. If it remains valid, the waves reflect off the ILR into 
long-trailing waves, which propagate out towards the CR. At the edge of the 
CR's forbidden region the waves morph into short-trailing waves, which prop- 
agate back towards the ILR. As they approach the ILR k is predicted to grow 
without limit. In reality the wave is absorbed as it approaches the ILR and its 
energy dissipates as heat. Thus the tight-winding approximation predicts that 
short-leading waves gradually wind up into short-trailing waves, which heat the 
disc in the vicinity of the ILR. The unwinding of leading waves and winding-up 
of trailing waves is similar to what differential rotation would do to material 
arms. 

Similarly, the dispersion relation implies that short leading waves will prop- 
agate inwards from the OLR to the outer edge of the forbidden region around 
the CR, where they will transfer to the long branch of the dispersion relation 
and move back out as long leading waves. If the tight-winding approximation 
remains valid as they approach the OLR, they will morph into long trailing 
waves that propagate back inwards towards the CR and then return as short 
trailing waves that eventually thermalise at the OLR. 

As waves transfer from leading to trailing form near a Lindblad resonance, 
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the waves morph from elongated ridges of overdensity to compact blobs of over- 
density. If Toomre's 

is small enough, self-gravity imparts a sharp inward impulse to these blobs, 
so their density begins to rise. Simultaneously, differential rotation is shearing 
them into trailing waves, which propagate away from the resonance. Conse- 
quently, the waves that reach the forbidden zone around the CR have larger 
amplitude than the waves that left this zone earlier. We say the waves have 
been swing amplified. Our formulae do not predict this behaviour because 
they rely on the tight-winding approximation, which is invalid at the crucial 
moment. 

The gain of the swing amplifier is a sensitive function of Q and the parameter 

X=^^, (1.101) 
m 

where fccrit is defined by equation (jl.98p . The smaller X is, the more invalid 
the tight-winding approximation, and the smaller Q is, the cooler the disc. The 
disc is stable to axisymmetric disturbances only if Q > iQ Swing amplification 
by a factor > 10 is possible for Q < 1.5 and X < 3. 

Since a key phase in the life-cycle of waves that we just described, the waves 
are unlikely to satisfy the tight-winding approximation, it is natural to a sk about 



solutio ns of the fundamental equation for normal modes, eq. (jl.79l) . iToomre 



( 198ll ) reported results obtained in this way by his student T. Zang. These 



showed that when the growth rate of a mode is not large, the mode looks like 
an interference pattern between leading and trailing waves that differ only in 
their amplitude, the trailing wave having the larger amplitude. The larger the 
mode's growth rate is, the more the trailing waves dominates. This finding is 
consistent with swing amplification taking place as disturbances morph from 
leading to trailing. Another finding was that t he modes are essentia lly confined 
to the region between the ILR and the CR. Read fc Evand ( IQQSl ) solved for 



the normal modes of discs with power-law gravitatinal potentials. They also 
found that the amplitudes of modes are largest between the ILR and the CR. 
Their models had "cut-out" discs, that is discs whose surface density tapered 
to negligible values at both very small and very large radii. The growth rate 
of a mode depended strongly on whether the ILR lay inside the inner cut-out, 
that is in a region of low density. In this case inward propagating short-trailing 
waves can reflect off the inner edge of the disc into leading waves, thus closing 
a feedback loop, rather than being absorbed at the ILR. Thus solutions to our 
mode equation (|1.79l) lend support to the qualitative understanding of spiral 
structure provided by the dispersion relation for tightly- wound waves. 

The bottom line is that stellar discs are responsive dynamical systems be- 
cause they support waves that can be amplified by self gravity as they move 
through the disc. The degree of amplification, and therefore the disc's respon- 
siveness, increases sharply as the velocity dispersion falls towards the critical 



*We can show this from equation l|1.98|l by setting m = 0. 
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value at which Q = 1 and the disc becomes unstable to radial fragmentation. 
Much of the energy carried by the waves is thermalised in the vicinity of a 
Lindblad resonance. Thus the waves heat the disc and render it less responsive. 



1.4.3 Spiral structure and normal modes 

iLin fc Shul i 1966h hypothesised that spiral structure is a manifestation of a 
mildly unstable normal model of the stellar disc: they envisaged the amplitude 
of this mode stabilising at a finite value as a result of energy dissipation in 
interstellar gas. They developed the theory of density waves in the expectation 
that the normal modes of discs could be understood in terms of waves trapped 
between barriers in the same way that we picture the modes of a laser as standing 
waves trapped between the laser's end mirrors. It's now clear that the very 
influential Lin~Shu paradigm is based on a misunderstanding of disc dynamics. 
Waves in a stellar disc are heavily damped already at the level of stellar dynamics 
because they heat the disc at the Lindblad resonances. 

From a certain perspective the failure of the Lin-Shu hypothesis is per- 
plexing: stellar dynamics is governed by the coupled Poisson and collisionless 
Boltzmann equations. These equations are time-translation invariant, so on 
group-theoretic grounds their linearised forms must have a complete set of so- 
lutions with time dependence e''*, with v possibly complex. Unless the system 
is completely stable, the evolution from any initial condition will be dominated 
by the most rapidly growing normal mode. Hence the observations must reflect 
such modes. 

The problem with this argument is that it assumes that any initial configura- 
tion can be represented by a superposition of normal modes. In other words, the 
normal modes are assumed to be complete. The solutions to our normal-mode 
equation ()1.79|) are not complete because in deriving it we have used defective 
logic: equation (jl.78p is obtained by dividing both sides of 

{n-n-Lj)SU = n-^S<^n (1.102) 

by n • r2 — oj. This operation is legitimate only if n • Jl — a; 7^ 0. If we want our 
normal modes to be complete, we have to include the case n • — cj = and 
replace equation (|1.78p by 

- "'^/°^^^ (5$„(J) + c„(J)(5(n -n-uj). (1.103) 
n- ll — uj 



In the simpler but closely analogous case of an electrostatic plasma, van KampenI 



was able to show that by considering such singular DFs Poisson's equa- 



tion can be satisfied for any real value of w. In this way we obtain a much 
richer set of solutions than can be obtained from equation (|1.78p . All these 
van Kampen modes are stable, and they prove to be complete, whereas the 
solutions we would obtain from equation (|1.78p are incomplete and as such do 
not form a basis for a discussion of stability. 
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1.4.4 Driving spiral structure 

By counting faint stars in the outer reaches of both our Galaxy and the An- 
dromeda nebula, M31, it has been shown that the outer p arts of galaxies are a 



mass of stellar streanis and full of fain t satel lite galaxies (jMcConnachie et al 



20091: iBelokurov et all , l2006t iBell et allboOSl ). From studies of the internal dy 



namics of satellite galaxies, we know that these systems are heavily dominated 
by dark matter, so we must anticipate that the dark-matter distribution that 
surrounds a galaxy like the Milky Way is lumpy. When a lump of dark mat- 
ter sweeps through pericentre, its tidal field will launch a wave into the host 
galaxy's disc, which we know to be a responsive system. The classic example 
of this process in M51, which has a satellite galaxy, NGC 5195, near the end of 
one of its exceptionally strong spiral arms. Few galaxies have such a luminous 
satellite so near to them, so grand-design spirals like that of M51 are not 
prevalent. Most galaxies will be responding simultaneously to more than one 
much weaker stimulus, with the result that their spirals are both weaker and 
rather chaotic. 

A majority of spiral galaxies have bars at their centres. The figures of bars 
are known to rotate quite rapidly in that the CR of the bar's pattern lies at 
a radius that is ^ 1.2 times the bar's length. The rotating gravitational field 
of the bar must perturb the disc, and from the discussion of ^1.4.21 we would 
expect the surrounding disc to show spiral structure at the pattern speed of 
the bar that extends from near the end of the bar to the OLR. However, both 
N-body simulations and observations show that the disc's principal respon se is 
at a lower pattern speed than that of the bar ( Sparke fc SellwoodL 1988[ ). so 



bar excites spiral structure that lies inside its CR. This finding is consistent 
with the tendency of the solutions to the normal-mode equation (|1.79p to have 
significant amplitudes only inside CR. Crucially it implies that the response to 
the bar rotates more slowly than the bar, so the relative phases of the features 
is constantly changing. 



1.5 Conclusion 

The key approximation of stellar dynamics is that stellar systems are collision- 
less, so the actual motion is well approximated by motion in a smooth potential. 
Orbits in smooth potentials are mostly quasiperiodic, and when they are not, it 
is possible to construct a nearby Hamiltonian in which the same initial condi- 
tions yield quasiperiodic motion. Therefore quasiperiodic motion is an excellent 
starting point for stellar dynamics. 

Quasiperiodic orbits display an elegant structure that is captured by angle- 
action coordinates: each orbit is a three-torus. The angle variables specify where 
on its torus a star is, and they evolve linearly in time. If the frequencies are 
incommensurable (as they nearly always should be) a star's probability density 
should be independent of angle variables, so the distribution function depends 
only on the actions. The actions provide a geometrical quantification of the 
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orbit. 

Normally the true Hamiltonian will differ from the approximate one that 
admits angle-action variables. The small difference h between the true Hamil- 
tonian and the approximate one can be important if it generates forces that act 
in the same way over extended periods of time. In the vicinity of resonances this 
may happen, and then h may change the dynamics qualitatively. The overall 
impact on the dynamics of a galaxy is nevertheless likely to be small if the reso- 
nance is isolated. When several resonances are simultaneously active, however, 
chaos can be generated, and stars may slowly diffuse through phase space. This 
process is likely to be important for the secular evolution of barred galaxies, but 
our understanding of it is currently inadequate. 

The potential of a real galaxy is always fluctuating, and fluctuations are 
of fundamental importance because they alone permit stars to exchange en- 
ergy. No matter what their physical origin, fluctuations will drive the system 
towards unattainable thermodynamic equilibrium, especially by enhancing core- 
halo structure. Two-body relaxation is mostly due to Poisson fluctuations in the 
number of stars in large volumes, and is an exceedingly slow process in galaxies. 
Hence in galactic dynamics we focus on fluctuations due to the motion of mas- 
sive bodies (giant molecular clouds, spiral arms dark-halo lumps, star-clusters 
and satellite systems). Fluctuations cause stars to diffuse through action space, 
and this diffusion is observed in the solar neighbourhood. The diffusion coefS- 
cients can be calculated from the temporal power spectrum of the fluctuations 
or empirically determined from observations of solar-neighbourhood stars. The- 
ory and observation are reasonably consistent, but there is plenty of scope for 
tightening constraints. 

Spiral structure is an important source of fluctuations. Its dominant effect is 
the creation of transient resonances, which by first trapping and then releasing 
stars cause them to move from inside the corotation circle outwards, and vice 
versa. The random velocities of stars are not increased by such churning, but 
stars with similar ages but different metallicities are mixed up. In addition to 
churning the disc around corotation, spiral structure moves angular momentum 
outwards, from ILR to OLR, in the process heating the disc in the vicinity of 
the Lindblad resonances, especially the ILR. 

Spiral structure is not fully understood because it is an inherently non-linear 
and global phenomenon. Lin-Shu-Kalnajs density-wave theory is restricted to 
the linear case and assumes tightly-would arms to make the physics essentially 
local. It predicts that tightly-wound leading waves propagate through a portion 
of the disc, unwinding as they go, so they inevitably violate the tight-winding 
approximation. As the waves pass from leading to trailing form, they are am- 
plified by a process that lies beyond linear theory, and eventually their energy is 
thermalised at a Lindblad resonance by an analogue of Landau damping. The 
inaccuracy inherent in using the tight- winding approximation can be eliminated 
by solving the exact equation for normal modes. Such solutions confirm the basic 
picture derived with the tight-winding approximation but reveals a preference of 
spiral structure to lie inside CR. Unfortunately, van Kampen's work on electro- 
static plasmas implies that the solutions of the normal-mode equation are not 
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complete, so they do not provide a secure basis for understanding the dynamics 
of discs. 

Self-gravitating discs are responsive systems because any disturbance is liable 

to excite leading waves, which may amplify significantly as they morph into 
trailing waves. Since the amplification becomes weaker as the disc heats, a 
pure stellar disc becomes less responsive as its velocity dispersion rises. Gas is 
an essential ingredient of a spiral galaxy because (i) it dissipates the energy of 
spiral waves, (ii) through star formation it constantly replenishes the population 
of stars with low velocity dispersion as spiral structure increases the velocity 
dispersion of older stars, and (iii) it makes any spiral gravitational potential 
observationally conspicuous by forming dust lanes and luminous blue stars near 
its troughs. 

Galaxies live in the noisy environments of their dark halos, into which clumps 
with various masses are continually falling. As they pass through pericentre 
such lumps may excite spiral structure. Any spiral structure will quickly heat 
the disc. If the disc is relatively cold and therefore responsive, much more 
energy will be converted into heat than was imparted by the exciting lump: 
by shifting angular momentum out through the disc, spiral structure makes 
gravitational energy available for random motions, and thus increases the disc's 
entropy. It is important never to lose sight of the fact that a disc of stars on 
nearly circular orbits is occupying only a tiny fraction of the phase space that 
is energetically accessible to it. Any random process will scatter its stars into a 
broader distribution in phase space, and thus make it a hotter, thicker disc. 
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